sf$FID<-row.names(sf)
sf<-merge(sf, s.nr.p, by.x="FID", by.y="IN_FID")
sf<-subset(sf, sf$NEAR_DIST==0)
sf<-as.data.frame(sf)
## Road quality.
## Join with precinct data.
ir<-read.dbf("C:/Dropbox (Nall Research Group)/petition/data/gislayers/mpos/PA_NHS_2012--NALambertConformalConic.dbf")$dbf
ir$FID<-as.numeric(row.names(ir))-1
ir<-as.data.frame(ir)
ir<-merge(ir, ir.nr.p, by.x="FID", by.y="IN_FID")
ir$IRI_VN[ir$IRI_VN==0]<-NA
## Bridge inventory
nb<-readShapePoints("C:/Dropbox (Nall Research Group)/petition/data/gislayers/mpos/nbi--joinedcounty--NALambertConformalConic")
nb$FID<-row.names(nb)
nb<-as.data.frame(nb)
nb<-merge(nb, n.nr.p, by.x="FID", by.y="IN_FID")
nb$ITEM67[nb$ITEM67=="*"]<-NA
nb$ITEM67<-as.numeric(as.character(nb$ITEM67))
## Now generate the summary statistics for each database, by NEAR_FID in the precinct file.
## SF summary
sf.ply<-ddply(sf, .(NEAR_FID), summarize, sf.stoplocs=length(unique(stop_id)),
sf.totalstops=sum(totalstops), sf.numagency=length(unique(agency)), sf.agencynames=paste(unique(agency), collapse="; "))
ir.ply<-ddply(ir, .(NEAR_FID.y), summarize, ir.totlength=sum(Shape_Leng), ir.meaniri.wt=weighted.mean(IRI_VN, Shape_Leng, na.rm=T), ir.meaniri.unwt=mean(IRI_VN))
nb.ply<-ddply(nb, .(NEAR_FID), summarize, nb.totbridges=length(unique(FID)), nb.mean.deficient=mean(ITEM67%in%c(3,2)), nb.mean.rating=mean(ITEM67))
## Assign MPO info, convex hull and county to Ansolabehere-Rodden data using the sp package over() function
names(ar.pts)
ar.pts.over.mpo<-over(ar.pts, mp.poly)
ar.pts.over.ch<-over(ar.pts, ch.poly)
ar.pts.over.co<-over(ar.pts, co.poly)
ar.pts2<-spCbind(ar.pts, ar.pts.over.mpo)
ar.pts2<-spCbind(ar.pts2, ar.pts.over.ch)
ar.pts2<-spCbind(ar.pts2, ar.pts.over.co)
## Now ready to start analysis of the
ar.pts3<-merge(ar.pts2, sf.ply, by.x="FID", by.y="NEAR_FID", all.x=T)
ar.pts3<-merge(ar.pts3, ir.ply, by.x="FID", by.y="NEAR_FID.y", all.x=T)
ar.pts3<-merge(ar.pts3, nb.ply, by.x="FID", by.y="NEAR_FID", all.x=T)
ard<-as.data.frame(ar.pts3)
## Merge in the data from Sanchez (2006) and Gerber and Gibson (2009)
write.csv(ard, "C:/Dropbox (Nall Research Group)/petition/data/gislayers/mpos/arwithtranspo.csv",row.names=F)
## mpo.analysis.R
library(apsrtable)
library(ggplot2)
library(biglm)
library(plyr)
## dissimilarity index code
seg.dissim<-function(g1pop, g2pop){
tmp<-na.omit(cbind(g1pop, g2pop))
g1pop<-tmp[,1]
g2pop<-tmp[,2]
g1prop<-sum(g1pop)/sum(g1pop+g2pop)
totpop<-sum(g1pop+g2pop)
popi<-g1pop+g2pop
g1propi<-g1pop/(g1pop+g2pop)
D<-sum(popi*abs(g1propi-g1prop)/(2*totpop*g1prop*(1-g1prop)))
return(D)
}
## Still need to bring in the MPO representation data: which ones are malapportioned, per Sanchez.
sa<-read.csv("C:/Dropbox/petition/data/mpos/gerber/mpoajpswithsanchez.csv", header=T, as.is=T)
ar<-read.csv("C:/Dropbox (Nall Research Group)/petition/data/gislayers/mpos/arwithtranspo.csv", header=T, as.is=T)
cr<-read.csv("C:/Dropbox (Nall Research Group)/petition/data/mpos/changednames.sa.csv", header=T, as.is=T)
mp<-read.csv("C:/Dropbox (Nall Research Group)/petition/data/mpos/gerber/mpoajpswithsanchez.csv", header=T, as.is=T)
## Update names in the
sa$new.mpo.sanchez<-sa$sanchez.mpo
for(i in 1:nrow(cr)){
sa$new.mpo.sanchez[sa$new.mpo.sanchez==cr$ggsaname[i]]<-cr$gisname[i]
}
sa$sanchez.mpo.old<-sa$sanchez.mpo
sa$sanchez.mpo<-sa$new.mpo.sanchez
ar<-merge(ar, sa, by.x="mp.MPONAME", by.y="new.mpo.sanchez", all.x=T)
ar$prop.obama<-ar$obama_shar
ar<-ddply(ar, .(mp.MPONAME), transform, dissim.obama=seg.dissim(NDV, NRV))
## International roughness index
## First, calculate by MPO voting rule bias
ar$sanchez.antiurb.bias<-ar$sanchez.urb.wt-ar$sanchez.urb.unw
ar$sanchez.antiurb.bias.cut4<-cut(ar$sanchez.antiurb.bias,4)
ar$sanchez.antiurb.bias.cut5<-cut(ar$sanchez.antiurb.bias,5)
ar$sanchez.urb.unw.cut4<-cut(ar$sanchez.urb.unw, 4)
ar$sanchez.urb.unw.cut3<-cut(ar$sanchez.urb.unw, 3)
ar$dissim.obama.cut3<-cut(ar$dissim.obama, 3)
ar$lnpopsqmi07<-log(ar$POP2007/(ar$Shape_Area/1000000))
ar$lnpopsqmi07[is.infinite(ar$lnpopsqmi07)]<-NA
ar$ln.ir.meaniri.wt<-log(ar$ir.meaniri.wt)
ar$sf.totalstops.pc<-ar$sf.totalstops/ar$POP2007
ar$sf.totalstops.sqkm<-ar$sf.totalstops/(ar$Shape_Area/1000000)
ar$sf.stoplocs.sqkm<-ar$sf.stoplocs/(ar$Shape_Area/1000000)
lm.1<-lm(ir.meaniri.wt~sanchez.antiurb.bias+sanchez.urb.wt, data=ar)
lm.2<-lm(ir.meaniri.wt~sanchez.antiurb.bias+sanchez.urb.wt+prop.obama, data=ar)
lm.3<-lm(ir.meaniri.wt~sanchez.antiurb.bias*lnpopsqmi07+sanchez.urb.wt+prop.obama, data=ar)
ggplot(data=ar, aes(x=lnpopsqmi07, y=ln.ir.meaniri.wt))+
stat_smooth(method=lm)+
facet_grid(.~sanchez.antiurb.bias.cut5)
## Primary exploratory finding: areas where cities have the lion's share of the MPO board have better urban roads
## And population density is more strongly (inversely) related to road quality in places with
p<-ggplot(data=subset(ar, !is.na(sanchez.urb.unw.cut4)), aes(x=lnpopsqmi07, y=ln.ir.meaniri.wt))+
stat_smooth(method=lm)+
facet_grid(sanchez.antiurb.bias.cut4~sanchez.urb.unw.cut4)+
geom_rug()
ggsave(p, filename = "C:/Dropbox (Nall Research Group)/petition/figs/iri.dens.sanchez.urb.unw.pdf")
## Related exploratory finding: high-density areas are more likely t
p<-ggplot(data=subset(ar, !is.na(dissim.obama.cut3)), aes(x=lnpopsqmi07, y=ln.ir.meaniri.wt))+
stat_smooth(method=lm)+
facet_grid(.~dissim.obama.cut3)+
geom_rug()
ggsave(p, filename = "C:/Dropbox (Nall Research Group)/petition/figs/iri.dissim.obama.pdf")
## NBI Bridge Deficiency Data
## MPO malapportionment data
p<-ggplot(data=subset(ar, !is.na(sanchez.urb.unw.cut4)), aes(x=lnpopsqmi07, y=nb.mean.deficient))+
stat_smooth(method=lm)+
facet_grid(sanchez.antiurb.bias.cut4~sanchez.urb.unw.cut4)+
geom_rug()
ggsave(p, filename = "C:/Dropbox (Nall Research Group)/petition/figs/nbi.dens.sanchez.urb.unw.pdf")
## Same story with bridge deficiency: in areas with more partisan geographic segregation
## Bridges in high-density areas are more likely to be deficient.
p<-ggplot(data=subset(ar, !is.na(dissim.obama.cut3)), aes(x=lnpopsqmi07, y=nb.mean.deficient))+
stat_smooth(method=lm)+
facet_grid(.~dissim.obama.cut3)+
geom_rug()
ggsave(p, filename = "C:/Dropbox (Nall Research Group)/petition/figs/nbi.dissim.obama.pdf")
##
## GTFS data: stop frequency by precinct.  Include agency fixed effects.
lm.out<-lm(sf.stoplocs.sqkm~prop.obama+sf.agencynames, data=ar)
p<-ggplot(data=subset(ar, !is.na(sanchez.urb.unw.cut4)), aes(x=lnpopsqmi07, y=log(sf.stoplocs.sqkm)))+
stat_smooth(method=lm)+
facet_grid(sanchez.antiurb.bias.cut4~sanchez.urb.unw.cut4)+
geom_rug()
ggsave(p, filename = "C:/Dropbox (Nall Research Group)/petition/figs/sf.stoplocs.sqkm.sanchez.urb.unw.pdf")
p<-ggplot(data=subset(ar, !is.na(sanchez.urb.unw.cut4) & sf.totalstops.pc<=5), aes(x=lnpopsqmi07, y=sf.totalstops.pc))+
stat_smooth(method=lm)+
facet_grid(sanchez.antiurb.bias.cut4~sanchez.urb.unw.cut4)+
geom_rug()
ggsave(p, filename = "C:/Dropbox (Nall Research Group)/petition/figs/sf.totalstops.pc.sanchez.urb.unw.pdf")
p<-ggplot(data=subset(ar, !is.na(sanchez.urb.unw.cut3) & sf.totalstops.pc<=5 & sf.totalstops.sqkm>0), aes(x=lnpopsqmi07, y=sf.totalstops.sqkm))+
stat_smooth(method=lm)+
facet_grid(.~sanchez.urb.unw.cut3)+
geom_rug()
ggsave(p, filename = "C:/Dropbox (Nall Research Group)/petition/figs/sf.totalstops.sqmi.sanchez.urb.unw.pdf")
p<-ggplot(data=subset(ar, !is.na(sanchez.urb.unw.cut3) & sf.totalstops.pc<=5 & sf.totalstops.sqkm>0), aes(x=lnpopsqmi07, y=sf.totalstops.sqkm))+
stat_smooth(method=lm)+
facet_grid(.~sanchez.urb.unw.cut3)+
geom_rug()
ggsave(p, filename = "C:/Dropbox (Nall Research Group)/petition/figs/sf.totalstops.sqmi.sanchez.urb.unw.pdf")
p<-ggplot(data=subset(ar, !is.na(dissim.obama.cut3)), aes(x=lnpopsqmi07, y=sf.totalstops.sqkm))+
stat_smooth(method=lm)+
facet_grid(.~dissim.obama.cut3)+
geom_rug()
ggsave(p, filename = "C:/Dropbox (Nall Research Group)/petition/figs/sf.totalstops.sqkm.dissim.obama.pdf")
#####################################################################################
## Limiting transit analysis to precincts in the "convex hull" of the transit service
#####################################################################################
lm.out<-lm(sf.stoplocs.sqkm~prop.obama+sf.agencynames, data=ar)
p<-ggplot(data=subset(ar, !is.na(sanchez.urb.unw.cut4)&!is.na(ch.NAME)), aes(x=lnpopsqmi07, y=log(sf.stoplocs.sqkm)))+
stat_smooth(method=lm)+
facet_grid(sanchez.antiurb.bias.cut4~sanchez.urb.unw.cut4)+
geom_rug()
ggsave(p, filename = "C:/Dropbox (Nall Research Group)/petition/figs/sf.stoplocs.sqkm.sanchez.urb.unw.ch.pdf")
p<-ggplot(data=subset(ar, !is.na(sanchez.urb.unw.cut4) & sf.totalstops.pc<=5&!is.na(ch.NAME)), aes(x=lnpopsqmi07, y=sf.totalstops.pc))+
stat_smooth(method=lm)+
facet_grid(sanchez.antiurb.bias.cut4~sanchez.urb.unw.cut4)+
geom_rug()
ggsave(p, filename = "C:/Dropbox (Nall Research Group)/petition/figs/sf.totalstops.pc.sanchez.urb.unw.ch.pdf")
p<-ggplot(data=subset(ar, !is.na(sanchez.urb.unw.cut3) & sf.totalstops.pc<=5 & sf.totalstops.sqkm>0&!is.na(ch.NAME)), aes(x=lnpopsqmi07, y=sf.totalstops.sqkm))+
stat_smooth(method=lm)+
facet_grid(.~sanchez.urb.unw.cut3)+
geom_rug()
ggsave(p, filename = "C:/Dropbox (Nall Research Group)/petition/figs/sf.totalstops.sqmi.sanchez.urb.unw.ch.pdf")
p<-ggplot(data=subset(ar, !is.na(sanchez.urb.unw.cut3) & sf.totalstops.pc<=5 & !is.na(ch.NAME)), aes(x=lnpopsqmi07, y=sf.totalstops.sqkm))+
stat_smooth(method=lm)+
facet_grid(.~sanchez.urb.unw.cut3)+
geom_rug()
ggsave(p, filename = "C:/Dropbox (Nall Research Group)/petition/figs/sf.totalstops.sqmi.sanchez.urb.unw.ch.pdf")
p<-ggplot(data=subset(ar, !is.na(dissim.obama.cut3)), aes(x=lnpopsqmi07, y=sf.totalstops.sqkm))+
stat_smooth(method=lm)+
facet_grid(.~dissim.obama.cut3)+
geom_rug()
ggsave(p, filename = "C:/Dropbox (Nall Research Group)/petition/figs/sf.totalstops.sqkm.dissim.obama.pdf")
## Data within convex hulls defined by transit
ar.ch<-ar[!is.na(ar$ch.NAME),]
names(ar.ch)
ar.ch$POINT_X
head(ar.ch$long)
names(ar.ch)
ar.ch$ar.long
ar.ch$ar.lon
ar.ch$ar.lat
unique(ar.ch$ch.NAME)
ar.ch$sanchez.mpo
unique(ar.ch$sanchez.mpo)
names(ar.ch)
install.packages(c("arm", "bayesm", "BH", "bit64", "blockTools", "BMA", "car", "chron", "class", "coin", "cwhmisc", "data.table", "dataframes2xls", "deldir", "DEoptimR", "dplyr", "e1071", "formatR", "Formula", "gdata", "GEOmap", "geosphere", "ggmap", "ggplot2", "ggtern", "glmnet", "goftest", "gplots", "gridExtra", "gtools", "highr", "Hmisc", "inline", "knitr", "lme4", "manipulate", "mapproj", "maps", "maptools", "markdown", "Matrix", "mgcv", "mime", "mvtnorm", "optmatch", "plotrix", "pscl", "psych", "R6", "raster", "Rcmdr", "RcmdrMisc", "Rcpp", "RcppArmadillo", "RcppEigen", "RCurl", "reldist", "reporttools", "reshape2", "rgdal", "rgl", "rJava", "RPMG", "Rwave", "sandwich", "sem", "sp", "spam", "spatstat", "spdep", "spgwr", "splancs", "stringi", "tidyr", "TSP", "vcd", "VGAM", "WriteXLS", "zoo"))
install.packages(c("arm", "bayesm", "BH", "bit64", "blockTools",
rm(list=ls())
q()
rm(list=ls()
)
q()
cw<-read.csv("C:/Dropbox (Nall Research Group)/dream/data/chettyetal2014/onlinedata4.dta")
ms<-read.dta("C:/Dropbox (Nall Research Group)/dream/data/chettyetal2014/onlinedata4.dta")
library(foreign)
ms<-read.dta("C:/Dropbox (Nall Research Group)/dream/data/chettyetal2014/onlinedata4.dta")
ms<-read.dta("C:/Dropbox (Nall Research Group)/dream/data/chettyetal2014/onlinedata4.dta")
names(ms)
ar<-read.csv("C:/Dropbox (Nall Research Group)/petition/data/gislayers/mpos/arwithtranspo.csv", header=T, as.is=T)
cw<-read.csv("C:/Dropbox (Nall Research Group)/dream/data/crosswalks/crosswalks.csv", header=T, as.is=T)
ms<-read.dta("C:/Dropbox (Nall Research Group)/dream/data/chettyetal2014/onlinedata4.dta")
ar<-read.csv("C:/Dropbox (Nall Research Group)/petition/data/gislayers/mpos/arwithtranspo.csv", header=T, as.is=T)
unique(ar$STCOFIPS)
unique(cw$msa2013)
ar$fips_12_13
library(foreign)
cw<-read.csv("C:/Dropbox (Nall Research Group)/dream/data/crosswalks/crosswalks.csv", header=T, as.is=T)
ms<-read.dta("C:/Dropbox (Nall Research Group)/dream/data/chettyetal2014/onlinedata4.dta")
ar<-read.csv("C:/Dropbox (Nall Research Group)/petition/data/gislayers/mpos/arwithtranspo.csv", header=T, as.is=T)
cw$msa2013
msa00<-read.dbf("C:/Dropbox (Nall Research Group)/dream/data/msa2000.dbf")
library(foreign)
msa00<-read.dbf("C:/Dropbox (Nall Research Group)/dream/data/chettyetal2014/msa2000.dbf")
names(msa00)
head(sort(msa00$CBSA_ID))
msa00<-read.dbf("C:/Dropbox (Nall Research Group)/dream/data/chettyetal2014/msa2000.dbf")
names(cms)
cms<-read.dta("C:/Dropbox (Nall Research Group)/dream/data/chettyetal2014/onlinedata4.dta")
names(cms)
cms$msa
cw<-read.csv("C:/Dropbox (Nall Research Group)/dream/data/crosswalks/crosswalks.csv", header=T, as.is=T)
names(cw)
head(sort(cw$fips))
head(sort(ar$STCOFIPS))
cw$msa2013
armsa13<-merge(ar, cw, by.x="STCOFIPS", by.y="fips")
head(armsa13)
armsa13<-NULL
cw2<-cw[,c("fips", "msa2013")]
armsa13<-merge(ar, cw2, by.x="STCOFIPS", by.y="fips")
nrow(armsa13)
cms$cbsatitle1
cms$cbsatitle2
head(cms)
msa.dissim<-ddply(armsa13, .(msa2013), transform, dissim.obama=seg.dissim(P2008_D, P2008_R))
library(plyr)
msa.dissim<-ddply(armsa13, .(msa2013), transform, dissim.obama=seg.dissim(P2008_D, P2008_R))
## Dissimilarity index code
seg.dissim<-function(grp1, grp2){
tmp<-cbind(grp1, grp2)
cat(dim(tmp), sum(grp1+grp2==0))
tmp<-tmp[tmp[,1]+tmp[,2]!=0,]
print("ran")
tmp<-na.omit(tmp)
print(nrow(tmp))
if(is.null(nrow(tmp))) return(NA)
if(nrow(tmp)==1) return(NA)
g1pop<-tmp[,1]
g2pop<-tmp[,2]
g1prop<-sum(g1pop)/sum(g1pop+g2pop)
totpop<-sum(g1pop+g2pop)
popi<-g1pop+g2pop
g1propi<-g1pop/(g1pop+g2pop)
D<-sum(popi*abs(g1propi-g1prop)/(2*totpop*g1prop*(1-g1prop)))
return(D)
}
msa.dissim<-ddply(armsa13, .(msa2013), transform, dissim.obama=seg.dissim(P2008_D, P2008_R))
msa.dissim
dim(msa.dissim)
names(msa.dissim)
tapply(msa.dissim$dissim.obama, msa.dissim$msa2013, unique)
armsa13$fips
names(cw)
names(cms)
head(cms)
## Read in CZ-MSA Crosswalk table
library(foreign)
library(plyr)
## Dissimilarity index code
seg.dissim<-function(grp1, grp2){
tmp<-cbind(grp1, grp2)
cat(dim(tmp), sum(grp1+grp2==0))
tmp<-tmp[tmp[,1]+tmp[,2]!=0,]
print("ran")
tmp<-na.omit(tmp)
print(nrow(tmp))
if(is.null(nrow(tmp))) return(NA)
if(nrow(tmp)==1) return(NA)
g1pop<-tmp[,1]
g2pop<-tmp[,2]
g1prop<-sum(g1pop)/sum(g1pop+g2pop)
totpop<-sum(g1pop+g2pop)
popi<-g1pop+g2pop
g1propi<-g1pop/(g1pop+g2pop)
D<-sum(popi*abs(g1propi-g1prop)/(2*totpop*g1prop*(1-g1prop)))
return(D)
}
cw<-read.csv("C:/Dropbox (Nall Research Group)/dream/data/crosswalks/crosswalks.csv", header=T, as.is=T)
cms<-read.dta("C:/Dropbox (Nall Research Group)/dream/data/chettyetal2014/onlinedata4.dta")
ar<-read.csv("C:/Dropbox (Nall Research Group)/petition/data/gislayers/mpos/arwithtranspo.csv", header=T, as.is=T)
## Merge MSA fips onto AnsRod data using the crosswalks.
cw2<-cw[,c("fips", "msa2013")]
armsa13<-merge(ar, cw2, by.x="STCOFIPS", by.y="fips")
msa.dissim<-ddply(armsa13, .(msa2013), transform, dissim.obama=seg.dissim(P2008_D, P2008_R))
names(msa.dissim)
tapply(msa.dissim$dissim.obama, msa.dissim$msa2013, unique)
cw<-read.csv("C:/Dropbox (Nall Research Group)/dream/data/chettyetal2014/List1.csv", header=T, as.is=T)
cms<-read.dta("C:/Dropbox (Nall Research Group)/dream/data/chettyetal2014/onlinedata4.dta")
ar<-read.csv("C:/Dropbox (Nall Research Group)/petition/data/gislayers/mpos/arwithtranspo.csv", header=T, as.is=T)
names(cw)
cw$stfips
cw$stcofips<-1000*cw$stfips+cw$cofips
cw$stcofips
armsa13<-merge(ar, cw, by.x="STCOFIPS", by.y="fips")
armsa13<-merge(ar, cw, by.x="STCOFIPS", by.y="stcofips")
names(armsa13)
msa.dissim<-ddply(armsa13, .(msa2013), transform, dissim.obama=seg.dissim(P2008_D, P2008_R))
names(armsa13)
armsa13$cbsafips
msa.dissim<-ddply(armsa13, .(cbsafips), transform, dissim.obama=seg.dissim(P2008_D, P2008_R))
names(msa.dissim)
tapply(msa.dissim$dissim.obama, msa.dissim$msa2013, unique)
tapply(msa.dissim$dissim.obama, msa.dissim$cbsatitle, unique)
names(cms)
names(cms)
head(cms)
sort(unique(cms$msa))
tapply(msa.dissim$dissim.obama, msa.dissim$cbsatitle, unique)
tapply(msa.dissim$dissim.obama, msa.dissim$cbsafips, unique)
tapply(msa.dissim$dissim.obama, as.character(msa.dissim$cbsafips), unique)
msa.dissim<-tapply(msa.dissim$dissim.obama, as.character(msa.dissim$cbsafips), unique)
msa.dissim
## Read in CZ-MSA Crosswalk table
library(foreign)
library(plyr)
## Dissimilarity index code
seg.dissim<-function(grp1, grp2){
tmp<-cbind(grp1, grp2)
cat(dim(tmp), sum(grp1+grp2==0))
tmp<-tmp[tmp[,1]+tmp[,2]!=0,]
print("ran")
tmp<-na.omit(tmp)
print(nrow(tmp))
if(is.null(nrow(tmp))) return(NA)
if(nrow(tmp)==1) return(NA)
g1pop<-tmp[,1]
g2pop<-tmp[,2]
g1prop<-sum(g1pop)/sum(g1pop+g2pop)
totpop<-sum(g1pop+g2pop)
popi<-g1pop+g2pop
g1propi<-g1pop/(g1pop+g2pop)
D<-sum(popi*abs(g1propi-g1prop)/(2*totpop*g1prop*(1-g1prop)))
return(D)
}
cw<-read.csv("C:/Dropbox (Nall Research Group)/dream/data/chettyetal2014/List1.csv", header=T, as.is=T)[1:1883,]
cms<-read.dta("C:/Dropbox (Nall Research Group)/dream/data/chettyetal2014/onlinedata4.dta")
ar<-read.csv("C:/Dropbox (Nall Research Group)/petition/data/gislayers/mpos/arwithtranspo.csv", header=T, as.is=T)
cw$stcofips<-1000*cw$stfips+cw$cofips
## Merge MSA fips onto AnsRod data using the crosswalks.
armsa13<-merge(ar, cw, by.x="STCOFIPS", by.y="stcofips")
msa.dissim<-ddply(armsa13, .(cbsafips), transform, dissim.obama=seg.dissim(P2008_D, P2008_R))
names(msa.dissim)
msa.dissim<-tapply(msa.dissim$dissim.obama, as.character(msa.dissim$cbsafips), unique)
msa.dissim<-data.frame(cbsafips=names(msa.dissim), dissim.obama=as.numeric(msa.dissim))
msa.dissim
cms<-merge(cms, msa.dissim, by.x="msa", by.y="cbsafips")
plot(cms)
names(cms)
plot(cms$s_rank_8082~cms$dissim.obama)
cor(cms$s_rank_8082, cms$dissim.obama, use="complete")
lines(lowess(cms$s_rank_8082~cms$dissim.obama))
names(cms)
pdf(file="C:/Dropbox/dream/figs/chettymsascores.obamadissim.pdf", paper="special", height=5, width=7)
plot(cms$e_rank_b~cms$dissim.obama, xlab="Obama Dissimilarity Index", ylab="Absolute Mobility")
lines(lowess(cms$e_rank_b~cms$dissim.obama))
dev.print(device=pdf, file="C:/Dropbox/dream/figs/chettymsascores.obamadissim.pdf")
pdf(file="C:/Dropbox/dream/figs/chettymsascores.obamadissim.pdf", paper="special", height=5, width=7)
plot(cms$e_rank_b~cms$dissim.obama, xlab="Obama Dissimilarity Index", ylab="Absolute Mobility")
lines(lowess(cms$e_rank_b~cms$dissim.obama))
dev.off()
plot(cms$s_rank_8082~cms$dissim.obama, xlab="Obama Dissimilarity Index", ylab="Chetty Rank-Rank Slope")
lines(lowess(cms$s_rank_8082~cms$dissim.obama))
plot(cms$e_rank_b~cms$dissim.obama, xlab="Obama Dissimilarity Index", ylab="Absolute Mobility")
lines(lowess(cms$e_rank_b~cms$dissim.obama))
cms$e_rank_b
lines(loess(cms$e_rank_b~cms$dissim.obama))
plot(cms$e_rank_b~cms$dissim.obama, xlab="Obama Dissimilarity Index", ylab="Absolute Mobility")
lines(lowess(cms$e_rank_b~cms$dissim.obama))
cms$e_rank_b
lines(lowess(cms$e_rank_b~cms$dissim.obama, na.rm=T))
lines(lowess(cms$e_rank_b~cms$dissim.obama, na.omit=T))
plot(cms$e_rank_b~cms$dissim.obama, xlab="Obama Dissimilarity Index", ylab="Absolute Mobility")
lines(lowess(cms$dissim.obama, cms$e_rank_b, na.rm=T))
?loess
plot(cms$e_rank_b~cms$dissim.obama, xlab="Obama Dissimilarity Index", ylab="Absolute Mobility")
lines(loess(cms$e_rank_b~cms$dissim.obama, na.action=na.omit))
View(cms)
plot(cms$e_rank_b~cms$dissim.obama, xlab="Obama Dissimilarity Index", ylab="Absolute Mobility")
lines(loess(cms$e_rank_b~cms$dissim.obama, na.action=na.exclude))
plot(cms$e_rank_b~cms$dissim.obama, xlab="Obama Dissimilarity Index", ylab="Absolute Mobility")
loess(cms$e_rank_b~cms$dissim.obama, na.action=na.exclude)
predict(loess(cms$e_rank_b~cms$dissim.obama, na.action=na.exclude))
lm.rank<-lm(cms$s_rank_8082~cms$dissim.obama)
summary(lm.rank)
lm.rank<-lm(cms$s_rank_8082~cms$dissim.obama, weights=cms$n_ige_rank_8082)
summary(lm.rank)
lm.abs<-lm(cms$e_rank_b~cms$dissim.obama, weights=cms$n_ige_rank_8082)
summary(lm.abs)
tmp<-tmp[tmp[,1]+tmp[,2]!=0,]
plot(cms$s_rank_8082~cms$dissim.obama, xlab="Obama Dissimilarity Index", ylab="Chetty Rank-Rank Slope")
abline(lm.rank)
predict(lm.rank)
pdf(file="C:/Dropbox (Nall Research Group)/dream/figs/chettymsascores.obamadissim.pdf", paper="special", height=5, width=7)
## Rank-based mobility measure (higher numbers indicates parent income more predictive)
lm.rank<-lm(cms$s_rank_8082~cms$dissim.obama, weights=cms$n_ige_rank_8082)
lm.abs<-lm(cms$e_rank_b~cms$dissim.obama, weights=cms$n_ige_rank_8082)
plot(cms$s_rank_8082~cms$dissim.obama, xlab="Obama Dissimilarity Index", ylab="Chetty Rank-Rank Slope")
abline(lm.rank)
## Absolute income mobility measure
plot(cms$e_rank_b~cms$dissim.obama, xlab="Obama Dissimilarity Index", ylab="Absolute Mobility")
dev.off()
dev.off()
dev.off()
cms$dissim.obama
pdf(file="C:/Dropbox (Nall Research Group)/dream/figs/chettymsascores.obamadissim.pdf", paper="special", height=5, width=7)
## Rank-based mobility measure (higher numbers indicates parent income more predictive)
lm.rank<-lm(cms$s_rank_8082~cms$dissim.obama, weights=cms$n_ige_rank_8082)
lm.abs<-lm(cms$e_rank_b~cms$dissim.obama, weights=cms$n_ige_rank_8082)
plot(cms$s_rank_8082~cms$dissim.obama, xlab="Obama Dissimilarity Index", ylab="Chetty Rank-Rank Slope")
abline(lm.rank)
## Absolute income mobility measure
plot(cms$e_rank_b~cms$dissim.obama, xlab="Obama Dissimilarity Index", ylab="Absolute Mobility")
dev.off()
pdf(file="C:/Dropbox (Nall Research Group)/dream/figs/chettymsascores.obamadissim.pdf", paper="special", height=5, width=7)
par(mfrow=c(1,2), las=1, cex=0.7, pty="s")
## Rank-based mobility measure (higher numbers indicates parent income more predictive)
lm.rank<-lm(cms$s_rank_8082~cms$dissim.obama, weights=cms$n_ige_rank_8082)
lm.abs<-lm(cms$e_rank_b~cms$dissim.obama, weights=cms$n_ige_rank_8082)
plot(cms$s_rank_8082~cms$dissim.obama, xlab="Obama Dissimilarity Index", ylab="Chetty Rank-Rank Slope")
abline(lm.rank)
## Absolute income mobility measure
plot(cms$e_rank_b~cms$dissim.obama, xlab="Obama Dissimilarity Index", ylab="Absolute Mobility")
abline(lm.abs)
dev.off()
cms$s_rank_8082
View(cms)
rm(list=ls())##clear workspace
install.packages(c("xlsx", "foreign", "stringr", "sandwich", "lmtest", "rms", "survey", "xtable", "lattice", "gplots", "graphics", "ggplot2", "car", "fossil", "maptools", "sp", "foreign", "zipcode", "reshape2", "plotrix"))
install.packages(c("xlsx", "foreign", "stringr", "sandwich",
lmtest", "rms", "survey", "xtable", "lattice", "gplots", "graphics", "ggplot2", "car", "fossil", "maptools", "sp", "foreign", "zipcode", "reshape2", "plotrix"))
install.packages(c("xlsx", "foreign", "stringr", "sandwich", "lmtest", "rms", "survey", "xtable", "lattice", "gplots", "graphics", "ggplot2", "car", "fossil", "maptools", "sp", "foreign", "zipcode", "reshape2", "plotrix"))
install.packages(c("xlsx", "foreign", "stringr", "sandwich",
install.packages(c("xlsx", "foreign", "stringr", "sandwich", "lmtest", "rms", "survey", "xtable", "lattice", "gplots", "graphics", "ggplot2", "car", "fossil", "maptools", "sp", "foreign", "zipcode", "reshape2", "plotrix"))
install.packages(c("xlsx", "foreign", "stringr", "sandwich", "lmtest", "rms", "survey", "xtable", "lattice", "gplots", "graphics", "ggplot2", "car", "fossil", "maptools", "sp", "foreign", "zipcode", "reshape2", "plotrix"))
install.packages(c("xlsx", "foreign", "stringr", "sandwich", "lmtest", "rms", "survey", "xtable", "lattice", "gplots", "graphics", "ggplot2", "car", "fossil", "maptools", "sp", "foreign", "zipcode", "reshape2", "plotrix"))
library(xlsx)
library(foreign)
library(stringr)
library(xlsx)
library(foreign)
library(stringr)
library(sandwich)
library(lmtest)
library(rms)
library(survey)
library(xtable)
library(lattice)
library(gplots)
library(graphics)
library(ggplot2)
library(car)
library(fossil)
library(maptools)
library(sp)
library(foreign)
library(zipcode)
library(reshape2)
library(plotrix)
as.n<-as.numeric
as.c<-as.character
as.f<-as.formula
len<-length
##Set your working directory here. Also store the file path to your working directory as the string object "wd".
wd<-"C:/Dropbox/geogaffect/replication package/final_upload_to_JOP"
setwd(wd)
##create file path for saving tables and figures
output<-paste(wd, "/output/", sep="")
### CONJOINT ANALYSIS
### Figure 2 in main text
### Figures 17-29 in Online Appendix
##Produce all results from conjoint analysis
source("conjoint.estimation.R")
### QUICKFIRE PAIRED COMPARISON ANALYSIS
### Figure 3 in main text
### Figures 3-16 in Online Appendix
source("quickfire.estimation.R")
### HEAT MAP ANALYSIS
### Figure 4 in main text
source("heatmap.estimation.R")
### NEIGHBORHOOD QUALITY PROXY ANALYSIS
### Display correlations between proxies of neighborhood quality and violent crime rates, county level discussed in Footnote 23
source("neighborhood.quality.proxy.analysis.R")
### FEASIBILITY OF SORTING ANALYSIS
### FIGURE 5 in main text
### FIGURES 32 - 55 in Online Appendix
source("feasibility.esimtation.R")
### OBSERVED MOVING BEHAVIOR OF SURVEY RESPONDENTS
### Figure 6 in main text
### Figure 56 in Online Appendix
source("move.behavior.estimation.R")
### COMMUNITY PREFERENCE FACTORIAL EXPERIMENT
### Figure 2 in Online Appendix
source("factorial.experiment.estimation.R")
##OTHER APPENDIX MATERIALS
##Tables in Online Appendix A
source("other.appendix.estimation.R")
update.packages()
y
q()
